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I .  INTRODUCTION 


This  investigation  is  part  of  an  ongoing  research  program  to  improve 
the  structural  response  modeling  capabilities  at  the  Terminal  Ballistics 
Division  of  the  Ballistic  Research  Laboratory.  Among  the  computer  programs 
used  to  model  structural  response  at  the  BRL,  the  ADINA  finite  element  code 
has  been  used  intensively  since  its  acquisition  in  1978.  It  is  a  well 
documented 1 »3  and  widely  used^»5,6  general  purpose  program  for 
calculating  the  static  and  dynamic  responses  of  complex  structures,  with  the 
capability  of  modeling  both  geometric  and  material  nonlinearities.  In 
particular,  concerning  the  elastic-plastic  behavior  of  materials,  the  code 
employs  the  Prandtl-Reuss  theory  of  plastic  flow^  with  either  the  von 
Mises  yield  condition  or  the  Drucker-Prager  yield  condition.  Using  the  von 
Mises  condition,  linear  isotropic  hardening  or  kinematic  hardening  can  be 
modeled,  while  only  elastic-perfectly  plastic  behavior  can  be  modeled  with 
the  Drucker-Prager  condition.  Although  linear  hardening  is  adequate  for 
situations  where  small  amounts  of  plastic  strain  are  involved,  its  use  in 
predicting  the  large  strain  response  of  materials  exhibiting  nonlinear 
plastic  behavior  can  result  in  considerable  error.  Hence,  in  order  to 
improve  the  modeling  of  strain  hardening  in  the  ADINA  code,  the  task  of 
incorporating  the  Besse ling-White  sublayer  model  into  the  ADINA  plasticity 
formulation  was  undertaken. 


1.  K.-J.  Bathe,  "ADINA  -  A  Finite  Element  Program  for  Automatic  Dynamic 
Incremental  Nonlinear  Analysis,"  Report  82448-1,  Acoustic  and  Vibration 
Lab,  MIT,  Dept  of  Mechanical  Engineering,  Sep  75  (rev.  Nov  79). 

2.  K.-J.  Bathe,  "Static  and  Dynamic  Geometric  and  Material  Nonlinear 
Analysis  Using  ADINA,"  Report  82448-2,  Acoustic  and  Vibration  Lab,  MIT, 
Dept  of  Mechanical  Engineering,  May  76  (rev.  May  77). 

3.  M.D.  Snyder  and  K.-J.  Bathe,  "Formulation  and  Numerical  Solution  of 
Thermo-Elastic-Plastic  and  Creep  Problems,"  Report  82448-3,  Acoustic  and 
Vibration  Lab,  MIT,  Dept  of  Mechanical  Engineering,  Jun  77. 

4.  K.-J.  Bathe  (ed.),  "Applications  Using  ADINA,"  Proceedings  of  the  ADINA 
Conference  August  1977,  Report  82448-6,  Acoustic  and  Vibration  Lab, 

MIT,  Dept  of  Mechanical  Engineering,  Aug  77. 

5.  K.-J.  Bathe  (ed.),  "Nonlinear  Finite  Element  Analysis  and  ADINA,"  Pro¬ 
ceedings  of  the  ADINA  Conference  August  1979  .  Report  82448-9,  Acoustic 
and  Vibration  Lab,  MIT,  Dept  of  Mechanical  Engineering,  Aug  79. 

6.  K.-J.  Bathe  (ed.),  "Nonlinear  Finite  Element  Analysis  and  ADINA,"  Pro- 
ceedings  of  the  3rd  ADINA  Conference,  MIT,  10-12  Jun  81  or  Computers 

&  Structures,  13,  No.  $-6,  1981. 

7.  R.  Hill,  The  Mathematical  Theory  of  Plasticity,  Oxford  University  Press, 
London,  1950,  pp.  38-45. 
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The  Besseling-White  model,  also  called  the  mechanical  sublayer  model, 
has  been  employed  for  some  time  in  a  number  of  structural  response 
codes**"!!  and  in  particular  in  the  PETROS  series  of  shell  response 
codes!^~!5  developed  for  the  BRL.  This  model  simulates  nonlinear 
strain  hardening  by  sunning  the  stresses  from  an  array  of  elastic- 
perfectly  plastic  elements,  called  sublayers,  for  the  material  stress  at  a 
point.  Nonlinear  hardening  is  modeled  in  a  piece-wise  linear  manner  by  each 
sublayer  yielding  at  a  different  stress  intensity,  with  the  sublayer  yield 
strength  being  determined  by  a  polygonal  fit  to  the  uniaxial  stress-strain 
curve. 


B.  H.A.  Balmer  and  E.A.  Witmer,  "Theoretical-Experimental  Correlation  of 
Large  Dynamic  and  Permanent  Deformations  of  Impulsively-Loaded  Simple 
Structures,"  Tech  Docu  Rpt  No.  FDL-TDR-64-108,  Jul  64,  Air  Force  Flight 
Dynamics  Lab,  Wright-Patterson  AFB,  Ohio. 

9.  J.M.  Santiago,  H.L.  Wisniewski  and  N.J.  Huffington,  Jr.,  "A  User's 
Manual  for  the  REPSIL  Code,"  BRL  Rpt  No.  1744,  Oct  74,  USA  Ballistic 
Res  Lab,  APG,  MD  (AD  A003176). 

10.  B.  Hunsaker,  Jr.,  D.K.  Vaughn  and  J.A.  Stricklin,  "A  Comparison  of  the 
Capability  of  Four  Hardening  Rules  to  Predict  a  Material's  Plastic 
Behavior,"  Second  National  Congress  on  Pressure  Vessels  and  Piping, 

ASHE,  San  Francisco,  CA,  23-27  Jun  75. 

11.  R.W.H.  Wu  and  E.A.  Witmer,  "Analytical  and  Experimental  Studies  of  Non¬ 
linear  Transient  Responses  of  Stiffened  Cylindrical  Panels,"  AIAA 
Journal,  Vol.  13,  No.  9,  Sep  75,  pp.  1171-1178. 

12.  L.  Morino,  J.W.  Leech  and  E.A.  Witmer,  "PETROS  2:  A  Finite-Difference 
Method  and  Program  for  the  Calculation  of  Large  Elastic-Plastic 
Dynamically-Induced  Deformation  of  Multilayer  Variable-Thickness  Shells" 

BRL  Contract  Rpt  No. 12,  Dec  69,  USA  Ballistic  Res  Lab,  APG,  MD  (AD  708774). 

13.  S.  Atluri,  E.A.  Witmer,  J.W.  Leech  and  L.  Morino,  "PETROS  3:  A  Finite- 
Difference  Method  and  Program  for  the  Calculation  of  Large  Elastic- 
Plastic  Dynamically-Induced  Deformations  of  Multilayer  Variable- 
Thickness  Shells,"  BRL  Contract  Rpt  No.  60,  Nov  71,  USA  Ballistic 

Res  Lab,  APG,  MD  (AD  890200L) . 

14.  S.D.  Pirotin,  B.A.  Berg  and  E.A.  Witmer,  "PETROS  3.5:  New  Developments 
and  Program  Manual  for  the  Finite-Difference  Calculation  of  Large 
Elastic-Plastic  Transient  Deformations  of  Multilayer  Variable-Thickness 
Shells,"  BRL  Contract  Rpt  No.  211,  Feb  75,  USA  Ballistic  Res  Lab, 

APG,  MD  (AD  A007215) . 

15.  S.D.  Pirotin,  B.A.  Berg  and  E.A.  Witmer,  "PETROS  4:  New  Developments 
and  Program  Manual  for  the  Finite-Difference  Calculation  of  Large 
Elastic-Plastic,  and/or  Viscoelastic  Transient  Deformations  of 
Multilayer  Variable-Thickness  (1)  Thin  Hard-Bonded,  (2)  Moderately-Thick 
Hard-Bonded,  or  (3)  Thin  Soft-Bonded  Shells,"  BRL  Contract  Rpt  No.  316, 

Sep  76,  USA  Ballistic  Res  Lab,  APG,  MD  (AD  B014253L). 


The  sublayer  model  has  been  programmed  in  the  ADINA  code  as  an 
additional  material  option  (model  12).  Presently,  it  is  available  for  use 
only  with  the  explicit  (central  difference)  time  integration  method.  In 
order  to  insure  that  the  sublayer  model  is  correctly  implemented  in  the  code, 
parallel  calculations  have  been  performed  using  the  sublayer  model,  the 
kinematic  hardening  model,  and  the  isotropic  hardening  model  to  represent  a 
linear  strain  hardening  material.  The  correlation  of  deflections  using  the 
sublayer  model  and  the  kinematic  hardening  model  is  excellent.  The 
correlation  between  the  sublayer  model  and  the  isotropic  hardening  model 
although  somewhat  poor,  as  is  to  be  expected,  is  still  reasonable. 

This  report  describes  the  formulation  of  the  sublayer  model  and  its 
implementation  in  the  ADINA  code.  It  also  gives  the  results  of  the 

comparison  between  the  sublayer  model,  the  ADINA  kinematic  hardening  model, 
and  the  ADINA  isotropic  hardening  model.  The  material  in  this  report  has 
been  published  in  preliminary  form  in  the  Proceedings  of  the  First  Chautauqua 
on  Finite  Element  Modeling. 


II.  FORMULATION  OF  THE  MECHANICAL  SUBLAYER  MODEL 


The  Besseling-White  or  mechanical  sublayer  model  originated  with  the 
concept  by  Duwez^>*®  of  modeling  nonlinear  plastic  hardening  and 

unloading  behavior  by  superposing  the  responses  from  a  continuous  aggregate 
of  elastic-perfectly  plastic  elements.  First,  White^  and  then 

Besseling20  in  a  more  general  setting  extended  Duwez's  one-dimensional 
model  to  three  dimensions  and  laid  the  basis  for  using  the  mechanical 
sublayer  model  as  a  computational  procedure  by  proposing  a  discrete  array  of 
elements  or  sublayers.  The  computational  use  of  the  mechanical  sublayer 

model  was  first  carried  out  at  MIT  by  Balmer  and  Witmer®  in  analyzing  the 
response  of  structures  using  the  finite  difference  method. 


16.  A.D.  Gupta,  J.M.  Santiago  and  H.L.  Wisniewski,  "An  Improved  Strain 
Hardening  Characterization  in  the  ADINA  Code  Using  the  Mechanical  Sub¬ 
layer  Concept,"  First  Chautauqua  on  Finite  Element  Modeling,  Harwichport, 
MA,  15-17  Sep  80,  pp.  335-351. 

17.  P.  Duwez,  "On  the  Plasticity  of  Crystals,"  Physical  Review,  Vol.  47, 

1935,  pp.  494-501. 

18.  H.F.  Bohnenblust  and  P.  Duwez,  "Some  Properties  of  a  Mechanical  Model 
of  Plasticity,"  J.  Appl.  Mech.,  Vol  15,  1948,  pp.  222-225. 

19.  G.N.  White,  Jr.,  "Application  of  the  Theory  of  Perfectly  Plastic  Solids 
to  Stress  Analysis  of  Strain  Hardening  Solids,"  Tech  Report  51,  Graduate 
Div  of  Applied  Math.,  Brown  University,  Aug  50. 

20.  J.F.  Besseling,  "A  Theory  of  Plastic  Flow  for  Anisotropic  Hardening  in 
Plastic  Deformation  of  an  Initially  Isotropic  Material,"  Report  S.  410, 
Nat.  Aero.  Rsch.  Inst.,  Amsterdam  (NLL),  1953. 
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The  mechanical  sublayer  model  assumes  that  the  physical  stress  a.  .  at  a 

point  in  the  material  can  be  decomposed  into  a  weighted  sum  of  N  sublayer 
a 

stresses  a. . 

ij 

N 


o. . 

ij 


£ 

a=l 


0) 


(1) 


with  u°  the  sublayer  weighting  factors.  Sublayers  are  assumed  to  be 

elastic-perfectly  plastic  and  to  experience  the  same  physical  strain  e..  . 

They  share  identical  elastic  constants,  but  have  distinct  values  for  their 
yield  strengths.  The  values  of  sublayer  yield  strengths  and  weighting 
factors  are  determined  from  a  polygonal  fit  to  the  uniaxial  stress-strain 
curve  (for  example,  see  Figure  2),  with  the  total  nun&er  of  sublayers  N  equal 
to  the  number  of  nonzero  slope  straight-line  segments  employed. 

We  can  briefly  describe  how  the  sublayer  model  simulates  nonlinear 
hardening  by  examining  the  case  of  uniaxial  stress.  Upon  initial  loading, 
all  the  sublayers  respond  elastically,  but  as  loading  proceeds  sublayers  will 
sequentially  become  plastic  as  they  each  reach  their  yield  limit.  Overall, 
this  results  in  progressively  diminishing  values  of  the  stress-strain  modulus 
as  the  total  stress  increases,  as  illustrated  in  Figure  1.  Upon  unloading 
and  reversal  of  loading,  all  the  sublayers  initially  become  elastic  and 
remain  elastic  until  the  sublayer  with  the  lowest  yield  stress  again  becomes 
plastic,  whereupon  the  process  repeats  itself  in  the  reverse  direction,  but 
with  one  significant  difference:  Straight-line  segments  are  now  twice  the 
length  of  the  corresponding  segments  on  the  initial  loading  curve  due  to  each 
sublayer  stress  having  to  change  from  a  given  value  of  the  yield  stress  to 
its  negative  value  rather  than  from  zero  to  the  given  value.  Thus,  the 
sublayer  model  exhibits  a  type  of  kinematic  hardening  behavior  first  proposed 
by  Masing21  that  has  subsequently  achieved  acceptance  in  modeling 
cyclic  loading. 22, 23 


21.  G.  Masing,  "Eigenspannungen  und  Verfestigung  be im  Messing,"  Proceedings 
2nd  Intern.  Congress  for  Appl.  Mech.,  Zurich,  Sep  1926. 

22.  Z.  Mroz,  "On  the  Description  of  Anisotropic  Workhardening,"  J.  Mech  Phys 
Solids,  Vol.  15,  1967,  pp.  163-175. 


23.  Z.  Mroz  and  N.C.  Lind,  "Simplified  Theories  of  Cyclic  Plasticity,"  Acta 
Mech.,  Vol.  22,  1975,  pp.  131-152. 


Figure  1.  Mechanical  sublayer  modeling  of  the  uniaxial  p  •  -«rve  upon 

initial  loading  and  reversal  of  loadi  >, 


A.  Derivation  of  tlie  Sublayer  Equations 


i 


i 


i 


t 


i 


An  incremental  plasticity  law  is  used  to  calculate  the  sublayer 
stresses.  For  a  given  sublayer  a,  the  sublayer  stress  increment  Aa? .  is 


related  to  the  elastic  strain  increment  Ae, 
law  for  the  material, 

A  0U.  = 


ij 


ij 

through  the  isotropic  Hooke's 


a  E  ®a  v  Aea  . 

1  +  v  |_  ij  1_2v  kk 


(2) 


where  F.  is  Young's  modulus  and  V  is  Poisson's  ratio.  The  strain  increment 

Ae.  .  is  the  sum  of  the  elastic  sublayer  strain  increment  A £..  and  the 

1J  Pa  1J 

plastic  sublayer  strain  increment  Ae. .  , 


ea  Pa 
Ae. .  =  Ae. .  +  Ae?. 
iJ  ij  i 3 


(3) 


where. as  already  mentioned,  the  total  strain  increment  Ae .  .  is  the  same  for 

all  sublayers  at  the  ooint .  Each  sublayer  is  assumed  to  obey  the  von  Mises 
yield  condition: 


where 


S?.  S?.  =  |  (aa) 

i]  ij  3  yJ 


a  is  the  sublayer  yield  stress  and 

y 


(4) 


ca  _  a  1  a  x 

^  •  “*  O .  .  “  o<  •  0.  . 
ij  ij  3  kk  ij 


(5) 


is  the  sublayer  deviatoric  stress.  Should  the  material  exhibit  strain  rate 
sensitivity,  the  yield  stress  for  each  sublayer  can  be  made  rate  dependent 
through  the  relation 


aa  =  aa  1  ♦  (e/Da)1/Na 

y  o  v 


(6) 


whereof  is  the  static  sublayer  yield  stress,  e2  is  the  second  invariant  of 

the  deviatoric  strain  rate  tensor,  and  0a,  N®  are  empirically  determined 

constants2*'  which  generally  are  assumed  to  be  identical  for  all  sub¬ 
layers.  To  complete  the  set  of  equations,  the  associative  flow  rule  is 
appended  in  order  to  determine  the  plastic  strain  increment: 


24.  N.  Perrone,  "Response  of  Pate-Sensitive  Frames  to  Impulsive  Load,"  J. 
Engr.  Mech.  Div.,  ARCE  Vol .  47,  No.  EMI,  Proceedings  Paper  No.  7890T~ 
Feb  71. 
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i??.  =  i\a  s?. 

1J  1J 


(7) 


where  AX  is  the  usual  plastic  flow  parameter  that  is  adjusted  to  keep  the 
stress  state  at  the  yield  value.  Combining  equations  (2),  (1),  and  (6)  we 
obtain  the  incremental  equation  used  by  the  sublayer  model  to  determine  the 
sublayer  stresses: 

-  ik  [%  *  i^k  4Ekk  s?j  <» 

where  for  convenience  we  replace  AX01  by 


AX 


'a 


(9) 


Equations  (1),  (4),  (5),  (8),  and,  when  treating  a  rate  sensitive 

material,  equation  (6)  are  the  basis  for  the  mechanical  sublayer  plasticity 
calculation.  When  the  sublayer  stress  rate  lies  inside  the  yield  surface,  as 

determined  by  (4),  AXa  in  (8)  is  set  equal  to  zero;  otherwise  AX°  is  deter> 

mined  by  the  requirement  that  the  resulting  stress  state  satisfy  the  yield 

condition.  The  details  of  determining  the  sublayer  weighting  factors  oia  and 

the  static  sublayer  yield  stresses  a°  from  the  polygonal  fit  to  the  uniaxial 

O 

stress^strain  curve  are  described  in  References  8  and  12. 

8 .  Computational  Procedure 

The  computational  procedure  assumes  that  at  a  given  time  step  or  cycle  n 
in  the  code's  computational  algorithm  the  strain  increments  Ae^  from  the 
previous  to  the  current  cycle  and  the  sublayer  stresses  a?,  (n-1)  at  the 
previous  cycle  are  known  for  all  Gaussian  integration  points.  Using  this 
information  the  sublayer  model  computes  first  the  current  sublayer  stresses 
a?.(n)  and  from  these  the  current  physical  stresses  (n)  bv  means  of  the 
algorithm^  which  we  now  outline. 


25.  J.M.  Santiago,  "Formulation  of  the  Large  T>ef lection  Shell  Equations  for 
Use  in  Finite^nif ference  Structural  Response  Codes,"  RRL  Rpt  No.  1571, 
Feb  72,  USA  Ballistic  Res  Lab,  APG,  MD  (AD  740742). 
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Starting  with  the  sublayer  having  the  smallest  value  of  the  yield  stress, 

a  trial  elastic  stress  increment  lB.  .  is  computed  bv  assuming  that  no  plastic 

ct 

flow  occurs  during  the  increment  so  that  AA  *0  and  equation  (A)  becomes: 


Jij  =  1+v  [Aeij  +  l-2v  Aekk  6ij  J  * 


Adding  the  elastic  stress  increment  to  the  previous  sublayer  stress,  a  trial 
sublayer  stress  for  the  current  cycle  is  obtained 

T  g 

°ij  =  °ij  Cn-1!>  +  A°ij  •  UD 

Writing  the  von  Mises  yield  function  for  the  sublayer  as 


_a  _a  1 


2  ,_ou 


<T(oC.)  =  a.,  o'*.  -  i  (a'r.)  -  ^  (an 
ij  i]  i]  3  v  kk'  3  y 


the  trial  sublayer  stress  is  checked  bv  substitution  in  the  von  Mises  yield. 

ct  Ta  ^ 

If  4>  (a..)  <0,  then  the  addition  of  the  elastic  stress  increment  Ao. . 

results  in  a  stress  state  inside  or  on  the  yield  surface,  and  hence  the  trial 

sublayer  stress  has  given  the  correct  sublayer  stress  for  the  current  cycle: 

°ij  <">  »  ^ “j  ■  03) 

T 

If,  however,  <p  ( a „)  >  O.then  yielding  occurs,  and  according  to  (8)  the 
trial  sublayer  stress  needs  to  be  corrected  to  obtain  the  current  sublayer 


stress: 


°ij(n)  =  °ij  "  ^ “  SyO1"1)  • 


To  determine  the  flow  parameter  AAa  ,  equation  (14)  is  inserted  in  the  von 
Mises  yield  function  (12)  and,  In  order  to  Insure  that  the  resulting  sublayer 
stress  state  is  on  the  yield  surface,  the  function  is  set  equal  to  zero;  l.e* 
^a(cr^j(n))  ■  0.  This  determines  a  quadratic  equation  in  AA**  (Reference 

25), which  is  solved  for  the  smallest  positive  root: 


where 


A  =  S?.(n-1)  S^Cn-l) 


(15) 


B  =  5°.  S“.(n-1) 

ij  ij 

c  -  *“(5“  )  . 


In  some  instances  when  the  strain  increment  is  too  large,  the  discrimi¬ 
nant  under  the  radical  in  equation  (15)  can  become  negative  due  to  the  trial 
sublayer  stress  causing  a  state  too  far  outside  the  yield  surface.  In  such 
instances  a  subincrement al  iterative  procedure  described  in  Reference  9  is 
employed  to  reduce  the  strain  increment  to  a  number  of  subincrementa  in  order 

to  obtain  a  real  solution  for  . 


Once  the  plastic  flow  parameter  AAa  is  determined,  the  sublayer  stress 

a.,  (n)  is  computed  from  (14)  and  is  stored  for  the  next  computational  cycle. 
1J 

This  orocedure  is  repeated  until  all  the  current  sublayer  stresses  at  the 
given  integration  point  are  found,  and  then  these  are  summed  with  the 
appropriate  weighting  factors  to  determine  the  physical  stress  at  that  point 
for  the  current  cycle 

N 

a.j(n)  =  a?j(n)  .  (16) 

o=l 


C.  Program  Implementation 


The  ADINA  code  is  structured  to  accept  nonlinear  material  models  not 
currently  included  in  the  material  library.  Ko  changes  to  the  program  are 
necessary  except  for  inserting  the  user  supplied  material  model  subroutine  in 
terms  of  the  appropriate  source  program  call  variables  described 
in  Reference  1. 
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The  sublayer  model  has  been  inserted  as  material  model  12  (subroutine 
EL3D12),  which  is  intended  for  triaxial  states  of  stress.  As  presently 
programmed  the  sublayer  model  can  only  be  used  for  explicit  time  integration. 
The  subroutine  calculates  the  current  stresses  from  the  strain  increments  and 
the  previous  stresses  as  explained  earlier,  but  has  not  as  yet  been 
programmed  to  calculate  the  stress-strain  matrix  (subroutine  MIDEP3)  needed 
to  determine  the  stiffness  matrix  required  for  implicit  time  integration. 
Because  the  sublayer  stresses  at  each  cycle  must  be  stored  in  order  to 
perform  the  stress  calculation  for  the  subsequent  cycle,  provision  for 
storage  must  be  made  in  accordance  with  the  number  of  sublayers  used.  This 
means  that  when  more  than  one  sublayer  is  employed  the  sublayer  model  will 
require  more  storage  than  the  ADINA  kinematic  hardening  model. 

III.  DYNAMIC  RESPONSE  ANALYSIS  OF  PLATE 


In  order  to  evaluate  the  ability  of  the  mechanical  sublayer  model  as 
implemented  in  the  ADINA  code  to  reproduce  the  dynamic  response  of  structures 
with  nonlinear  material  characteristics,  calculations  have  been  undertaken  to 
compare  code  predictions  with  known  experimental  results.  Some  preliminary 
results  from  this  investigation  are  now  available  concerning  a  comparison 
between  the  sublayer  model  and  the  existing  ADINA  hardening  models. 
Employing  each  plasticity  model,  the  code  was  used  to  solve  a  problem 
involving  the  large  transient  deflection  of  an  impulsively  loaded, 
clamped-edge,  circular  flat  plate. 

The  problem  is  taken  from  an  experimental  report^**  and  concerns  a 
plate  6.35mm  thick  with  a  radius  of  63.5mm,  made  of  2024-0  aluminum.  This 
material  exhibits  a  considerable  degree  of  nonlinear  hardening  which  can  be 
closely  approximated  by  the  sublayer  model  as  illustrated  in  Figure  2. 


However,  for  purposes  of  comparing  the  sublayer  model  (model  12)  to  the 
isotropic  hardening  model  (model  8)  and  the  kinematic  hardening  model  (model 
9),  we  need  only  consider  the  following  material  responses:  purely  elastic, 
elastic-perfectly  plastic,  and  elastic-linear  strain  hardening.  Hence,  for 
this  purpose  the  following  values  of  the  material  properties  were  used: 


0  (yield  stress 

y 

E  (Young's  modulus) 

E<j>  (strain  hardening  modulus) 
V  (Poisson’s  ratio) 
p  (mass  density) 


*  85.5  MPa 

•  73.7  GPa 
“  6.85  GPa 
»  1/3 

-  2775  Kg/m3 


The  deformation  was  initiated  by  subjecting  the  plate  to  a  uniform  impulse 
velocity  of  53.09  m/s. 


26.  E.A.  Witraer,  F.  Merlis  and  S.D.  Pirotin,  "Experimental  Studies  of 

Explosively-Induced  Large  Deformations  of  Flat  Circular  2024-0  Aluminum 
Plates  with  Clamped  Edges  and  of  Free  Thin  Cylindrical  6061-T6  Shells," 
BRL  Contract  Rpt  No.  13&,  Jan  74,  USA  Ballistic  Res  Lab,  APG,  MD 
(AD  917518L). 
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Figure  2.  Tensile  uniaxial  stress-strain  data  for  2024-0  aluminum 
taken  from  Reference  26  apnroximated  by  a  number  of 
material  response  representations. 


STRESS  (ksi) 


By  passing  a  symmetry  plane  through  a  diameter,  only  half  the  plate  had 
to  be  treated.  Twenty  3-dimensional  brick  elements  were  used  to  model  half 
the  plate,  as  illustrated  in  Figure  3.  Each  element  used  16  nodes  except  for 
those  bounding  the  pole  which  used  only  12  nodes,  giving  a  total  of  142  nodes 
for  the  problem.  Two  Gaussian  integration  points  were  employed  along  each 
coordinate  direction,  making  a  total  of  eight  points  per  element.  All  nodes 
along  the  circular  boundary  were  assumed  to  be  fixed.  Explicit  time 
integration  was  used  for  all  calculations  and  the  total  Lagrangian 
formulation  was  used  for  the  nonlinear  geometric  analysis. 

A.  Elastic  Comparison 

The  purely  elastic  response  was  modeled  by  increasing  the  value  of  the 

yield  stress  a  to  a  level  insuring  elastic  behavior  only. (See  Figure  2.)  The 

elastic  calculations  were  performed  using  first  a  linear  and  then  a  nonlinear 
geometric  analysis. 

For  the  linear  analysis,  the  elastic  response  of  the  sublayer  model 
(model  12)  was  compared  with  the  response  of  the  ADINA  linear  analysis  model 
(model  1)  and  the  elastic  response  of  the  ADINA  strain  hardening  model 
(either  model  8  or  9,  since  the  isotropic  and  kinematic  models  coincide  for 
elastic  response).  The  numerical  results  from  the  three  models  were  found  to 
coincide  exactly.  This  is  illustrated  in  Figure  4  for  the  history  of  the 
deflection  at  the  pole. 

For  the  geometric  nonlinear  analysis  the  elastic  responses  of  the 
sublayer  model  and  the  strain  hardening  model  were  compared.  Again,  it  was 
found  that  the  models  gave  coincident  numerical  results,  as  shown  In  Figure  4. 

Hence  for  both  linear  and  nonlinear  geometric  analyses,  the  sublayer 
model  and  the  ADINA  models  are  elastically  consistent.  Moreover,  it  is 
observed  that  for  the  given  problem,  nonlinear  geometric  effects  are  already 
substantial  at  levels  of  maximum  deflection  in  the  order  of  60Z  of  the  plate 
thickness. 


B.  Blaatic-Perfectly  Plastic  Comparison 

For  the  elastic-perfect ly  plastic  comparison,  the  yield  stress  was 

returned  to  the  value  of  0^  *  83.5  MPa  and  the  strain  hardening  modulus  E-p 

was  set  equal  to  zero,  as  shown  in  Figure  2.  The  sublayer  model  (with  only 
one  sublayer),  the  isotropic  hardening  model,  and  the  kinematic  hardening 
model  were  compared  using  the  nonlinear  geometric  analysis. 


Excellent  correlation  between  the  results  of  the  three  models  was  found 
although  they  no  longer  coincided  numerically.  As  portrayed  in  Figure  4, 
the  deflection  history  curves  for  the  models  overlap.  Hence,  for  elastic- 
perfectly  plastic  behavior,  the  sublayer  mode 1  and  the  ADINA  plasticity  models 
give  the  same  results  physically.  Also,  comparing  the  elastic-perfectly 
plastic  response  with  the  purely  elastic  response,  we  see  the  dissipative 
effects  of  plastic  flow  in  the  reduced  and  smoothed  deflection  peaks  at  late 
times. 


18 


Figure  3.  Finite  element  modeling  of  half  the  plate  showing  a  typical  element 

and  the  direction  of  polar  deflection. 


C.  Elastic-Linear  Strain  Hardening  Compariaon 

To  model  strain  hardening  behavior,  the  strain  hardening  modulus  was 
reset  at  the  value  of  ■  6.85  GPa.  Two  sublayers  were  used  with  the 

yield  stress  of  the  second  sublayer  raised  to  a  value  high  enough  to  insure 
its  elastic  behavior,  so  that  only  the  first  sublayer  could  become  plastic. 
This  insured  that  the  sublayer  model  would  exhibit  linear  hardening  behavior, 
so  that  a  meaningful  comparison  could  be  made  with  the  standard  ADINA  models. 
The  nonlinear  geometric  analysis  option  was  used  for  the  comparison. 

The  correlation  between  deflections  from  the  sublayer  model  and  the 
kinematic  hardening  model  is  still  very  good,  as  illustrated  in  Figure  5,  but 
the  curves  no  longer  overlap  as  they  did  for  the  perfectly  plastic  case.  The 
agreeawnt  between  the  values  of  the  stress  components  at  Gaussian  integration 
points  ranged  in  the  order  of  one  to  three  significant  figures. 

On  the  other  hand,  the  deflection  using  the  isotropic  hardening  model 
was  only  in  fair  agreement  with  the  deflections  from  the  other  two  models, 
although  still  physically  acceptable.  However,  this  is  hardly  surprising 
since  the  isotropic  model  bases  its  plasticity  calculation  on  the  concept  of 
an  expanding  yield  surface  rather  than  a  shifting  yield  surface  as  do  the 
other  models.  Hence,  whenever  a  significant  amount  of  plastic  flow  reversal 
occurs,  the  models  can  be  expected  to  predict  different  responses. 

Calculations  were  also  performed  using  four  sublayers  to  model  the 
stress-strain  curve  as  shown  in  Figure  2.  However,  for  the  level  of  loading 
imposed  on  the  plate,  most  of  the  plastic  flow  was  experienced  by  the  first 
sublayer,  so  that  the  results  were  graphically  indistinguishable  from  those 
of  the  two-sublayer  model. 


IV .  CONCLUSIONS 

The  mechanical  sublayer  model  has  been  iaq>lemented  in  the  ADINA  finite 
element  structural  response  code,  extending  its  capability  to  model 
elastoplaatic  constitutive  behavior  to  nonlinear  kinematically  hardening 
materials.  The  model  is  formulated  to  treat  triaxial  states  of  stress.  It 
is  currently  programmed  for  use  in  both  the  linear  and  nonlinear  geometric 
analyses  of  transient  response  employing  explicit  tism  integration.  The 
sublayer  model  subroutine  has  been  exercised,  and  calculations  using  the 
sublayer  model  and  the  standard  ADINA  hardening  models  show  that  results  are 
in  close  agreement  and  appear  physically  reasonable. 

A  sublayer  model  for  biaxial  states  of  stress  is  currently  being  devised 
to  treat  plane  stress  problems.  When  ready,  this  model  will  also  be 
implemented  in  the  code.  It  now  appears  that  the  sublayer  model  can  be 
easily  reformulated  to  permit  implicit  time  integration.  Further  effort  is 
underway  to  validate  the  model  with  available  experimental  results  and  with 
predictions  from  reliable  finite  difference  response  codes,  such  as  REPSIL^ 
and  PETROS.**  Eventually,  the  sublayer  model  can  be  expected  to 

replace  the  existing  kinematic  hardening  model,  thereby  considerably 
extending  the  material  modeling  capabilities  of  the  ADINA  code. 


%  ■■ 
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Figure  4.  Comparison  of  the  deflection  histories  at  the  Dole  for  purely  elastic 
responses  and  an  elastic-perfectly  plastic  response. 


Figure  5.  Comparison  of  the  deflection  histories  at  the  pole  for  an  elastic-linear  hardening 
response  using  the  isotropic  hardening,  kinematic  hardening,  and  sublayer  models. 
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1.  Element  Group  Control  Card  (2014) 


Page  X.5  (Ref.  1) 


# 


Columns 

Variable 

Model 

57-60 

NPAR( 15) 

NCON 

65-68 

NPAR( 17) 

IDW 

77-80 

NPAR(20) 

Material 

Property 

Data  Cards 

a.  Material  number  card 

No 

change 

b.  Material  Property  Card 

Columns 

Variable 

1-10 

PROP ( 1 ,  N  ) 

11-20 

PR0P(2,N) 

3.  Sublayer  Data 

a. 

where 

NSUBL(ILAY) 

NUMMAT  -  NPAR(16) 

b.  (SIG,DZ(ILAY,ISB) 
ISG.D7. 

NSBL 
ILAY 

c.  (EPS, DZ( ILAY, ISB) 
EPS ,DZ 

NSBL 


12,  for  sublayer  model 
2,  for  sublayer  model 
(Number  of  sublayers  x  6)  +  7 

Page  X.ll  (Ref.  1) 

Page  X.ll  (Ref.  1) 

Young's  modulus,  F 
Poisson's  ratio,  v 


(NSDBLdLAY),  ILAY  -  1, NUMMAT) 


(1615) 


the  number  of  sublayers  in  layer  ILAY. 

Number  of  different  sets  of  material  properties. 

ISD  -  l.NSBL)  (5E15.6) 

*  the  magnitude  of  the  ordinate  of  the  1-dimensional  stress- 
strain  curve. 

■  Number  of  sublayers  in  the  layer  (material)  in  question. 

■  Material  in  question. 

ISB-1,  NSBL)  (5F.15.6) 

the  magnitude  of  the  abscissa  of  the  1-diraensional  stress- 
strain  curve. 

Number  of  sublayers  in  the  layer  (material)  in  question. 


-w 
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no  o  on  o  o  n  n  n  n  n  n  n  o  o  o  o  n  n  o  n  o  o  n  o 


EVALUATE  STRESS  INCREMENTS  AND  STRESSES 
MOOEL*12 

THE  FOLLOWING  VARIABLES  ARP  USED  IN  THIS  SUBROUTINE 

SIG  PPEVIOUS  STRESSES 
EPS  PREVIOUS  STRAINS 

STRESS  CURRENT  STRESSES  < TO  BE  CALCULATED) 

STRAIN  CURRENT  STRAINS  (GIVEN) 

IPEL  *  1*  MATERIAL  FLASTIC  (INITIAL  VALUE) 

*  ?.  MATERIAL  PLASTIC 
OELEPS  INCREMENTAL  STRAINS 

DELSIG  INCREMENTAL  STRESSES.  CALCULATED  ON  THE  ASSUMPTION 

OF  ELASTIC  BEHAVIOR  DURING  STRAIN  INCREMENT  (DELfPS) 
L  NO.  OF  SUR-INCRFMENTS 
LC  CHECK  ON  SUB-INCREMENTS 
PROP (l )  YOUNGS  MODULUS 
PROP (2 )  POISSONS  RATIO 
TR  TRIAL  STRESSES 

TC  CORRECTOR  STRFSSES 

TM  STRESS  PEP/SURLAYER 

DIMENSION  PROP  ( 1 )  lEPSm  .SIb(NSBL*M 

COMMON  /EL/  INO.ICOUNT.NPAR (20) .NUMEG.NEGl .NE GNL . IMASS. IDAMP . 

i  istat.ndof.nl  IN. ieig.imassn.idampn 

COMMON  /VAR/  NG.MOnFX.lUPDT.KSTEP.ITEMAA.IE0REF.ITE.KPRI, 

I  ipef.ifouit.ipri.kplotn.kplote 

COMMON  /MTM03D/  0 ( 36 ) »STPE SS (6 ) .STR A IN (ft ) • IPT .NEL . IPS 
COMMON  /MATER/  COEFF (6) .SIGMA (6).EfctT(6) 

COMMON  /MAlpPI/  (VSURL(A),SIGI0/(A.S>.ERS1DZ(*,5) 

COMMON  /DISOR/  DISQ (9) 

COMMON  /ELSTP/  TIME . IDTRF 

DIMENSION  DFLSlG(ft)  . DF.LFPS  (6)  .TR  (ft)  »TC  (ft)  ,TM  (6)  .STATE  (?) 
EQUIVALENCE  (NPARO)  .INDNL  ) 

DATA  STATE  /?H  E,ZH*P/ 

IPELD* IPEL 

IF ( IPT  .NE.  1 ) GOTO  30 

.  CALCULATION  OF  MATERIAL  CONSTANTS  . . . . 

YM«PROP ( 1 ) 

PV«PROP(?) 

CALL  MATPRO  (MATP.PROP) 

SI0M*SIGMA(1) 

FACT0R«1. 

AI*YM/(I,.PV) 

C1*A1/?. 

Al«AI/(l,-?.*PV) 

B1«A1*PV 

A1*A1-B1 

l.  calculate  incremental  strains 

30  00  32  1*1.6 
STRESS ( I ) *0. 0 
3?  CONTINUE 
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IFUNONL  .EO.  3 ) GOTO  36 

TOTAL  LAGRANGIAN 

00  34  1*1*6 

34  DFLEPSU)*STPAIN(I)-EPSm 
GOTO  40 

UPDATED  LAGPANGI AN 

36  00  36  1*1*6 
38  DELEPS (I ) *STPA IN ( I ) 

40  IPLAST*0 

00  250  ISt>*l  *NSBL 

STRESS  PEP  SUBLAYER  SI6(1SB*I> 

00  43  1*1*6 
TM(I) *SIG(ISH.I I 
43  CONTINUE 

KPRI*0  SKIP  CALCULATION  JUST  PRINT  STRESS 
IF  (KPRI  .£0.  0 ) GOTO  228 
L*  l 
f  L  *  1  • 

SIGPSGxSIGMA ( ISP)*SIGNA ( I Sb ) *F ACTOR*F ACTOR 

2.  CALCULATE  TRIAL  ELASTIC  STRESS  INCREMENTS 

PELSIG(1)*A1*DELEPS(1)  ♦  Hl*(DFLERS(2)*0ELEPSt3>> 

DELSIG (2) *A l*CELEPS (2)  ♦  PI* (DELfRS ( 1 ) ♦UELEPS(3> ) 
DELSIG(3)*Al*DFLfPS(3)  ♦  Bl* (DFLEPS (1) ♦DELEPS (2) > 
DELSIb(4)«Cl*DELFPS(4) 

DELSIG(5)*C1*DELEPS(5) 

DELSIG(b)*Cl*DELtPS(6) 

45  LC*1 

50  SUB  INC* 1 • /EL 

3.  CALCULATE  TRIAL  ELASTIC  STRESSES 
00  60  1*1*6 

TR<I)*TM< I > ♦DELSIG <I)*SUtINC 
60  CONTINUE 

TAUR*TR(1) ♦TR { 2) ♦TR13) 

4.  CHECK  THROUGH  VON  MISES-HENCKY  YIELO  CRITERIA 
YIELD  PHI 

C:Z*TP(1)**2>TP(2)**?^TR(3>**2*2.0*(TR(4)«*24TRC5)**?^TR(6)**2) 
1 (TADR**242.0*S1GMSO)/3.0 

TEST  YIELD  CONDITION 

HLAROAbO.O 

IF(CZ.LE.0.  .AND.  L  .Ea.UGOTO  215 
IF (CZ  *LE.  0 • ) GOTO  110 

5.  COMPUTE  NECESSARY  CORRECTOR 

SM*(TM(l)^TM(?)tTM(3))/3*0 
T  C ( 1 ) *TM ( 1 ) -SM 
TC ( 2 ) *TM (?) -SM 


o n  oon  noon  no  n  noon  o  o  o  o  o  o  o  o  o  o  non 


TC(3)*TM(3)-SM 
TC 141 *TM (4) 

TC (5) *TM (5) 

TC  (6 1  *TM (A) 

TAUC-TC  UWTC(2)*TC<3> 

6.  CALCULATE  PLASTICITY  PARAMETER  FROM  HUFF IN6T0N  MOOEL 

AZsTC (1)#*?*TC<?)#*2*TC(3)*#2*2*0*ITC 14) **2*TC (51 **2*TC (fcl *#?) 
1  -(TAUC*«?/3.0) 

P2*TCCl)*TRU)*TC(2l*TH(2>*TC(3i«TH(3)*2.0«(TPU)*TC(4>* 

1  T«(5)*TC  (S)  ♦TRU»)*TC(AI  )-TAUC*TAUR/3.0 

DlSCR*BZ*P2-A7*CZ 

TEST  A7 

IF  A2  IS  NEGATIVE  -  PRINT  ERROR  MESSAGE 
IF  A Z  IS  ZERO  -  SUP-INCREMENT 

IF  AZ  IS  POSITIVE  -  CONTINUE 

IF(AZ)  80  » 1  SO  *  1 00 
80  WRITE <P»R0> 

90  FORMAT  (»3*«4X«*AZ  NEGATIVE  AT  •/) 

GOTO  180 

TEST  DISCRIMINANT 

IF  DISCR  IS  NEGATIVE  SUH-lNCREMENT 
100  IF (DISCR  .LT.  0.01GCTO  150 
TEST  PZ 

IF  B2  IS  NEC-ATIVE  OR  7ERO  SUB-INCREMENT 

IF(BZ  .LE.  0.0) GOTO  150 
compute:  hlambda 
HLAMOA«CZ/(RZ«SORT(OISCRI  I 

7,  EVALUATE  SUBLAYER  RELIEVED  ELASTIC  STRESSES 

110  00  120  1*1*4 

TM ( I ) *TR ( I ) -ML AMDA*TC ( I ) 

120  CONTINUE 


CHECK  THE  SUB-INCREMENT  NUMREP 

IF(LC.EO.L)  GOTO  210 

LC*LC»l 

GO  TO  50 

MAKF  SUB-INCREMFNTS  SMALLER 


150  L*L*1 

EL*L 


CHECK  MAXIMUM  number  OF  ALLOWABLE  SUB-INCREMENTS 


uu  uuu  uu  u 


IFa  .Lt.  100) GOTO  45 
WRITE (6.170) 

170  FORMAT (*3STPESS  CALCULATION  UNSATISFACTORY (100  SUB- INCREMENTS) 
1  /) 

160  WRITE  16*3000)  ISB.L.LC.NEL. IPT.TIME 
3000  F0RMaT(*  SUBLAYER**.^,*  SUB-INCPEMMENT  L*».I4,»  LC**.I4. 

1*  ELEMENT  NUMBER** • 14* •  IPOINT** • U. •  TIME** .E15.8) 

WRITE (6*3005)  CZ.A7.bZ 

3005  FORMAT ( •  CZ**«E15.H«*  AZ**.E15.B.*  HZ**.E15.H) 

WRITE (6*3010)  (SIG(ISB.I) .1*1.6) 

3010  FORMAT  (  •  SIC.(ISB.I)  *«6Flb.H) 

WRITE (6.3020)  (TP ( I ) .1*1.4) 

3020  FORMAT ( •  TR  *«6E15.ft) 

WRITE (6.3021 )  (TC ( 1 ) . 1*1 .6) 

3021  F OPMAT ( •  TC  *«6E15.6) 

WRITE (6.3022)  (TM ( I ) .1*1.6) 

3022  FORMAT ( •  TM  *.6E15.0) 

WRITE (6.3023)  (OELSIG ( I ) • I *1 .6 ) 

3023  FORMAT ( •  OELSIG  *.6E15.6) 

STOP  • STRESS  calculation  UNSATISFACTORY* 

C  REACHED  PLASTIC  SOLUTION 

210  iplast.iplastm 

GOTO  225 

C  ELASTIC  TM  EQUALS  TRIAL  TR  PER  SUBLAYER 

215  00  220  1*1.6 
TM ( I ) *TH  !  I ) 

220  CONTINUE 

STRESS  ROTATION  IS  APPLIED  IN  LARGE  DISPLACEMENT /STRAIN 

225  IF ( INUNL.NE *3) GOTO  226 
0MEGAl«(>ISn(4)-0ISU<6) 

0MEGA2*DISD(5)-0ISD(M) 

OME6A3«OISD (T)-OISD(V) 

TM(l)*TM(l).siG(ISB.4)*0ME GA1*SIG( ISB.5) •0ME6A2 
TM(2)«TM(2)-S1G(ISB.4)»0ME6A1.SIG(ISB.6)*0MEGA3 
TM(3)*TM(3)-SIG(ISb.5)*0MEGA2-SIG(ISB.6)*0MEGA3 
TM(4)«TM(4) .,5*(0MEGA1*(SIG(ISB«2)-SIG(ISB.1) )♦ 

1  OMEGA 3*SIG ( ISB.5) ♦0MEGA2*SIG ( ISB.6) ) 

TM(5)*TM(5)..5*(0MEGA2*(SIG(IS6.3)-SIG(ISB.1) )♦ 

1  OMEGA  I *S IG ( ISB.6) -OMEGA 3*S 16 ( I SB. 4) ) 

TM(6)*TM(6) ♦  ,5* (OMEGA 3* (SIG(ISB«3)-SIG(ISB.2) ) — 

1  0MEGA2*S1G(ISR.4)-0MEGA1*SIG(1SB.5) ) 

226  IF(IUPOT  .NE.  0)GOTd  226 

STORE  SUBLAYER  STRESS  TMU)  IN  SIG(ISB.I)  FOR  NEXT  STER 

00  227  1*1.6 
SIG(ISfi.I)*TM(I) 

227  CONTINUE 

8.  CALCULATE  TOTAL  STRESSES  FROM  SUBLAYER  STRESS 

228  DO  230  1*1.6 

STRESS ( I ) -STRESS ( I ) *TM ( I ) *COEFF ( ISH ) 

230  CONTINUE 
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non  n  n  r>  o  n  <■>  n  r>  n  no 


2S0  CONTINUE 


.  CALCULATION  OF  ELASTIC-PLASTIC  STRESSES  . .  N  0 

IF (KPRI  ,E0.  0 ) GOTO  700 
IPELD»1 

IFdPLAST  ,6T.  OJIPELP*? 

UPDATING 

IF  < IUPUT.NF.O)GOTO  Gib 
I PFL* IPELD 
1*0  610  1*1.6 
610  f PS < I ) »ST«AIN ( I) 

615  RETURN 

PRINTING  OF  STPESSES 

700  IF  (IPPI.F0.0  .AND.  IPT.ht.l)  WRITE  (6.2100)  Nf L 

9.  CALCULATE  HYDROSTATIC  AND  OEVIATORlC  STRESSES 

SM* (STRESS ( 1 ) .STRESS (2) .STRESS (3) >/3.U 

SXsSTRtSS(l)-SR 

SY=STPFSS(?)-SR 

S?*STRESS<3)-SR 

F  TA*  .5  *  (SX«*2.SY**?.SZ**2>  ♦ 

1  STrE SS ( 4 ) **2  ♦  STRESS (b) **?  ♦  STMESS(6)**H 

FTA«SOPT(3.*FTA) 

IF  (INDNL.NE.?)  GO  TO  600 

IN  TOTAL  LAGPANGIAN  FOnHULaTICN  CALCULATE  CAUCHY  STRESSES 
IF  (IPRI.EO.O) 

1  WRITE  (6.2200)  IPT. STATE (IPELD) .STRESS. FTA.L.LC 
CALL  CAUCH3 

600  IF  (IPRl.EQ.O) 

1  WRITE  (6.2200)  IPT.STATF (IPELD) .STRESS.FTA.L.LC 

RETURN 

2100  EORHAT  (/I6H  ELEMENT  STRESS. VOX , 1 OHtQUI VaLENT/ 

1  131R  NUM/IPT  STATE  STHESS-XX  STRFSS-YY 

2S-Z2  STPFSS-XY  STRESS-XZ  STrESS-YZ  STRESS 

3L  LC  / I 4 ) 

2200  FORMAT  (h> . 1 2 . 1 X . A2 .6HL AST IC » I X .6E 1 4. 6. U.f 14 ,6.3X. 1 3.SX. 13) 


LIST  OP  SVMROLR 


a 

N  empirically  determined  constants  in  the  strain  rate  equation. 
Young's  modulus  (GPa). 
tangent  strain  hardening  modulus  (GPa). 
deviator ic  stress  tensor  at  the  ath  sublayer. 

Kronecker  delta. 

second  invariant  of  the  deviator ic  strain  rate  tensor, 
physical  strain  at  a  point, 
hydrostatic  comnonents  of  strain. 

Poisson's  ratio. 

mass  density  (Ke/m^). 

physical  stress  at  a  point  (MPa). 

hydrostatic  components  of  the  sublayer  stress  (MPa). 

(n)  sublayer  stress  at  the  ath  sublayer  at  the  current  cycle  (MPa). 
(n-1)  sublayer  stress  at  the  ath  sublayer  at  the  previous  cycle  (MPa), 
trial  sublayer  stress  (MPa), 
sublayer  yield  stress  (MPa), 
static  sublayer  yield  stress  (MPa), 
von  Mises  yield  function, 
sublayer  weighting  factor, 
strain  increment, 
plastic  flow  parameter, 
sublayer  stress  increment  (MPa), 
elastic  stress  increment  (MPa). 
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